Short-time dynamics and critical behavior of three-dimensional site-diluted Ising 

model 
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Monte Carlo simulations of the short-time dynamic behavior are reported for three-dimensional 
weakly site-diluted Ising model with spin concentrations p = 0.95 and 0.8 at criticality. In con- 
trast to studies of the critical behavior of the pure systems by the short-time dynamics method, 
our investigations of site-diluted Ising model have revealed three stages of the dynamic evolution 
characterizing a crossover phenomenon from the critical behavior typical for the pure systems to 
behavior determined by the influence of disorder. The static and dynamic critical exponents are 
determined with the use of the corrections to scaling for systems starting separately from ordered 
and disordered initial states. The obtained values of the exponents demonstrate a universal behavior 
of weakly site-diluted Ising model in the critical region. The values of the exponents are compared 
to results of numerical simulations which have been obtained in various works and, also, with results 
of the renormalization-group description of this model. 

PACS numbers: 64.60. De, 64.60.Ht, 61.43.-j, 75.40.Mg 



I. INTRODUCTION 

The investigation of critical behavior of disordered sys- 
tems remains one of the main problems in condensed- 
matter physics and excites a great interest, because all 
real sohds contain structural defects [l|, 0|- The struc- 
tural disorder breaks the translational symmetry of the 
crystal and thus greatly complicates the theoretical de- 
scription of the material. The influence of disorder is 
particularly important near the critical point where be- 
havior of a system is characterized by anomalous large 
response on any even weak perturbation. The descrip- 
tion of such systems requires the development of special 
analytical and numerical methods. 

The effects produced by weak quenched disorder on 
critical phenomena have been studied for many years (3|- 
[S]. According to the Harris criterion Q, the disorder 
affects the critical behavior only if a, the specific heat 
exponent of the pure system, is positive. In this case 
a new universal critical behavior, with new critical ex- 
ponents, is established. In contrast, when a < 0, the 
disorder appears to be irrelevant for the critical behav- 
ior. Only systems whose effective Hamiltonian near the 
critical point is isomorphic to the Ising model satisfy this 
criterion. 

A large number of publications is devoted to the study 
of the critical behavior of diluted Ising-like magnets by 
the renormalization-group (RG) methods, the numerical 
Monte Carlo methods, and experimentally (for a review 
see Refs. @, M SS; ^^^ [ii|). The ideas about replica 
symmetry breaking in the systems with quenched dis- 
order were presented in Refs. [l^, U^- A refined RG 
analysis of the problem has shown the stability of the 
critical behavior of weakly disordered three-dimensional 



systems with respect to the replica symmetry breaking 
effects [ij]. All obtained results confirm the existence 
of a new universal class of the critical behavior, which 
is formed by diluted Ising-like systems. However, it re- 
mains unclear, whether the asymptotic values of critical 
exponents are independent of the rate of dilution of the 
system, how the crossover effects change these values, 
and whether two or more regimes of the critical behavior 
exist for weakly and strongly disordered systems. These 
questions are the subjects of heated discussions @, Ha]. 

For critical dynamic systems, traditionally it is be- 
lieved that universal scaling behavior is realized in the 
long-time regime of dynamic evolution. However, at first 
in the paper [lg|, it was shown that systems, starting 
from macroscopic non-equilibrium initial states, demon- 
strate a universal scaling behavior on the macroscopic 
short-time stages of their dynamic process which is char- 
acterized by initial slip exponents and 9' for the re- 
sponse functions G{r, t, t') and the order parameter m{t) 
(magnetization for ferromagnetic systems) 



G{r,t,t')^{t/t'f, m{t)^t'> 



(1) 
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A remarkable property of this relaxation process is the 
increase of magnetization m{t) from a non-zero initial 
magnetization toq <C 1 at short times t < tcr ^ 
TOq z^'^)^ The initial rise of magnetization is 

changed to the well known decay m{t) ~ f-0/^'^ for 
t ^ tcr UM- The critical exponents 9 and 9' depend on 
the dynamic universality class [17| and have been calcu- 
lated by the RG method for a number of dynamic mod- 
els [18| such as the model with a non-conserved order 
parameter [1^, [3 (model A) , the model with an order 
parameter coupled to a conserved density [20] (model C), 
and the models with reversible mode coupling [2ll | (mod- 
els E, F, G, and J). The universal scaling behavior of the 
initial stage of the critical relaxation for pure systems 
has been verified by extensive numerical simulations |22| - 
|25| . The developed method in these papers of short-time 



critical dynamics gives the possibility to determine both 
the static critical exponents v, /3, and dynamic critical 
exponents z and 9' in the macroscopic short-time regime 
of the critical relaxation. However, a number of publi- 
cations devoted to the numerical study of disorder influ- 
ence on non-equilibrium critical relaxation by the short- 
time dynamic method is surprisingly little. We know 
the papers |26iT2^] in which the non-equilibrium criti- 
cal dynamics of the three-dimensional (3D) site-diluted 
Ising ferromagnets with quenched point-like defects is in- 
vestigated and the values of the initial slip exponent 0' 
27 [ and an exponent Ca for autocorrelation function [2^ 
are determined for systems with different spin concen- 
trations. The obtained universal value 9' = 0.10(2) is in 
good agreement, as it is insisted in \2IT\, with the RG es- 
timate for 9' = 0.0868 calculated in the two-loop approx- 
imation in Ref. |29| with the use of e-expansion method, 
where £ = 4 — d, with d is the spatial dimension. How- 
ever, some assumptions introduced during investigations 
in Ref. [23] and discussed below precludes from consent- 
ing to this value 9' = 0.10(2) as confirmation of the RG 
estimate validity. In our paper [31| the integrated Monte 
Carlo simulations of the short-time dynamic behavior are 
reported for 3D Ising and XY models with long-range cor- 
related disorder at criticality, in the case corresponding 
to linear defects. Both static and dynamic critical expo- 
nents are determined for systems starting separately from 
ordered and disordered initial states. The obtained val- 
ues of the exponents are in good agreement with results 
of the field-theoretic description of the critical behavior 
of these models in the two-loop approximation [30| • 

In the present paper, we numerically investigate the 
short-time critical dynamics with a non-conserved order 
parameter (model A) [131 in the 3D site-diluted Ising sys- 
tems with spin concentrations p — 0.95 and 0.8. In the 
following section, we introduce the 3D Ising model with 
quenched point-like defects and scaling relations for the 
short-time critical dynamics. In Sec. HI, we derive the 
critical short-time dynamics in Ising systems starting sep- 
arately from ordered and disordered initial states. Crit- 
ical exponents obtained under these two conditions with 
the use of the corrections to scaling are compared. The 
final section contains analysis of the main results, their 
comparison to results of other investigations and our con- 
clusions. 



Si fixed at the lattice sites and assuming values of ±1. 
Nonmagnetic impurity atoms form empty sites. In this 
case, occupation numbers pi assume the value or 1 and 
are described by the distribution function 



P{pi) = {l-p)5{p,)+p5{l-pi), 



(3) 



with p = 1 — c, where c is the concentration of the impu- 
rity atoms. 

In this paper we have investigated systems with the 
spin concentrations p = 0.95 and 0.8. We have considered 
the cubic lattices with linear size L = 128. The Metropo- 
lis algorithm has been used in simulations. We consider 
only the dynamic evolution of systems described by the 
model A in the classification of Hohenberg and Halperin 
\Vt\ . The Metropolis Monte Carlo scheme of simulation 
with the dynamics of a single-spin fiips reflects the dy- 
namics of model A and enables us to compare obtained 
critical exponents 9' and z to the results of RG descrip- 
tion of the non-equilibrium relaxation of this model. 

According to the argument of Janssen, et al. [la| ob- 
tained with the RG method and e-expansion, one may 
expect a generalized scaling relation for the fcth moment 
the magnetization 

m^''^t,T,L,mo) ^ b-'"^/''m^'''> (b-H, b^/^r, b'^L, &^°mo) 

(4) 
is realized after a time scale imic which is large enough in 
a microscopic sense but still very small in a macroscopic 
sense. In Eq. (|3]), 6 is a spatial rescaling factor, /? and v 
are the well-known static critical exponents, and z is the 
dynamic exponent, while the new independent exponent 
xq is the scaling dimension of the initial magnetization 
Too and t = {T — Tc)/Tc is the reduced temperature. 

Since the system is in the early stage of the evolution 
the correlation length is still small and finite size prob- 
lems are nearly absent. Therefore we generally consider 
L large enough (L = 128) and skip this argument. We 
now choose the scaling factor b = t^/^ so that the main 
t-dependence on the right is cancelled. Applying the scal- 
ing form (d]) for fc = 1 to the small quantity t^°^^mo, one 
obtains 

m(t,r,mo) - mo/ F{t^^'''T,t'''>^'mo) (5) 

= TOo/(l + ati/'^V) + 0(T2,m2), 



II. DESCRIPTION OF THE MODEL AND 
METHODS 

We have considered the following 3D site-diluted fer- 
romagnetic Ising model Hamiltonian defined in a cubic 
lattice of linear size L with periodic boundary conditions: 



H = -J 



2^ PiPj 



JiJ' 



J' 



(2) 



where the sum is extended to the nearest neighbors, J > 
is the short-range exchange interaction between spins 



where 9 — (xq — P/i')/z has been introduced. It was 
shown in Ref. [ig that the critical exponents 9 and 
9' in Eq. ^ are related by the scaling relation 9' = 
9 + {2 — z — T])l z, therefore independent exponent is one 
of them {9 or 9'). For r = and small enough t and 
TO,o, the scaling dependence for magnetization ([5]) takes 
the form m{t) ^ t^ . The time scale of a critical initial 
increase of the magnetization is tcr ~ to,q ^'^°. However, 
in the limit of toq -^ 0, the time scale goes to infinity. 
Hence the initial condition can leave its trace even in the 
long-time regime. For illustration we give in Fig. 1 the 
time evolution of the magnetization m(t) from the initial 




10000 t(MCS/s) 



FIG. 1: Time evolution of tiie magnetization m{t) from the 
initial state with magnetization mo = 0.03 at Tc — 3.49948 
as a result of Monte Carlo simulation of samples with spin 
concentration p = 0.8 and with linear size L = 128. 



state with magnetization mo = 0.03 at Tc — 3.49948 as 
a result of Monte Carlo simulation of samples with spin 
concentration p = 0.8 and with linear size L — 128. 

If T ^ 0, the power law behavior is modified by the 
scaling function F{t^/'"^) with corrections to the sim- 
ple power law, which will be dependent on the sign of 
T. Therefore, simulation of the system for tempera- 
tures near the critical point allows to obtain the time- 
dependent magnetization with nonperfect power behav- 
ior and the critical temperature T^ can be determined by 
interpolation. 

For the site-diluted Ising model, we measured the time 
evolution of the magnetization determined as follows: 



m[t) 






(6) 



where angle brackets denote the statistical averaging, the 
square brackets are for averaging over the different impu- 
rity configurations, and Ns — pL^ is a number of spins in 
the lattice. Two other interesting observables in short- 
time dynamics are the second moment of magnetization 
m(2)(f). 



.(2) 



it) 






(7) 



and the auto-correlation function 



A{t) 



^J2pAit)S,{0) 



(8) 



As the spatial correlation length in the beginning of the 
time evolution is small, for a finite system of dimension 



d with lattice size L the second moment m^'^\t, L) ^ L'^. 
Combining this with the result of the scaling form in 
Eq. (j4]) for T = and b ~ t^/^, one obtains 

m(2)(t) ^ t-2/3/-V2)(l,t-V-L) ^t^^ (9) 



C2 



d-2- 



Furthermore, careful scaling analysis shows that the 
auto-correlation also decays with a power law |32| 



A(t)-r 



Ca = I 

Z 



(10) 



Thus, the investigation of the short-time evolution of sys- 
tem from a high-temperature initial state with toq — 
allows to determine the dynamic exponent z, the ratio of 
static exponents 13 /v, and the initial slip exponent 9' . 

Until now, a completely disordered initial state has 
been considered as starting point, i.e., a state of very 
high temperature. The question arises how a completely 
ordered initial state evolves, when heated up suddenly to 
the critical temperature. In the scaling form (|3]), one can 
skip besides i, also the argument rriQ — 1, 



,(fc) 



(i, r) = b-'^^'^m'^''^ (b-H, b^/^A . (11) 



The system is simulated numerically by starting with a 
completely ordered state, whose evaluation is measured 
at or near the critical temperature. The quantities mea- 
sured are m{t) and m^'^^t). With b = t^^^, one avoids 
the main t dependence in m^'^' (t) and for k = 1 one has 

m{t,T) = t-P^'''m{l,f^'''T) (12) 

= i-/^/"^ (l + at^/''''T + 0{t'^)] . 

For r = 0, the magnetization decays by a power law 
7n{t) ^ t^*^/"^. If T 7^ 0, the power law behavior is 
modified by the scaling function m{l, I^/'^^t). From this 
fact, the critical temperature Tc and the critical exponent 
P/vz can be determined. 

The scaling form of magnetization in Eq. (|12p is pre- 
sented as follows: 

In m(t,r) = (-/3/z/z) Ini + lnm(l, t^/^^V) (13) 

after differentiation with respect to r gives the power law 
of time dependence for the logarithmic derivative of the 
magnetization in the following form: 



dr\nm{t,T)l^^^t^''' 



(14) 



which allows to determine the ratio \/vz. On the basis of 
the magnetization and its second moment, the cumulant 



U2{t) 



.(2) 



(to) 2 



1 - t'^'^ 



(15) 



is defined. From its slope, one can directly measure the 
dynamic exponent z. Consequently, from an investiga- 
tion of the system relaxation from ordered initial state 
with mo = 1, the dynamic exponent z and the static ex- 
ponents /3 and v can be determined and their values can 
be compared to results of simulation of system behavior 
from disordered initial state with niQ ~ 0. 



III. MEASUREMENTS OF THE CRITICAL 

EXPONENTS FOR 3D SITE-DILUTED ISING 

MODEL 

We have performed simulations on three-dimensional 
cubic lattices with linear size L — 128, starting either 
from an ordered state or from a high-temperature state 
with zero or small initial magnetization. We would like to 
mention that measurements starting from a completely 
ordered state with the spins oriented in the same direc- 
tion (ttiq = 1) are more favorable, since they are much 
less affected by fluctuations, because the quantities mea- 
sured are rather big in contrast to those from a random 
start with mo = 0. Therefore, for careful determination 
of the critical exponents for 3D Ising model with spin 
concentrations p = 0.95 and 0.8, we begin to investigate 
the relaxation of this model from a completely ordered 
initial state. 



A. Evolution from an ordered state with mo = 1 

Initial configurations for systems with the spin con- 
centrations p = 0.95; 0.8 and with randomly distributed 
quenched pointlike defects were generated numerically. 
Starting from those initial configurations, the system 
was updated with Metropolis algorithm at the criti- 
cal temperatures T^ = 4.26267(4) for p — 0.95 and 
Tc = 3.49948(18) forp = 0.8, which have been deter- 
mined in our paper [33| using Monte Carlo simulation of 
the 3D site-diluted Ising model with different spin con- 
centrations and particularly with p — 0.95 and 0.8 in 
equilibrium state. At present investigation, simulations 
have been performed up to t = 1000 Monte Carlo steps 
per spin (MCS/s). We measured the time evolution of 
the magnetization m(i) and the second moment mS^>(t\ 
which also allow to calculate the time-dependent cumu- 
lant U2{t) inEq. ([T51). 

In Fig. 2 the magnetization m,(i) is plotted on a log-log 
scale for samples with spin concentrations p — 0.95 [Fig. 
2, a] and -p = 0.8 [Fig. 2, b] at T = Tc(p) (curves 1) and 
at T = Tc(j>) T AT with AT = 0.005 (curves 2 and 3). 
In Fig. 3, the logarithmic derivative of the magnetization 
9r lnm(i, t)|^^q with respect to r and in Fig. 4 the cu- 
mulant f/2(0 ^"^^ plotted on a log- log scale at T = Tc{p) 
for samples with spin concentrations p = 0.95 (a) and 
p = 0.8 (b), accordingly. The 9,- lnm(i,T)|^^p have 
been obtained from a quadratic interpolation between 
the three curves of time evolution of the magnetization in 



Fig. 2 for the temperatures T = Tc(p), T = Tc(p) T AT 
and taken at the critical temperature T^ — 4.26267(4) 
for samples with p = 0.95 and T^ = 3.49948(18) for 
p — 0.8. The resulting curves in Figs. 2-4 have been 
obtained by averaging over 6000 and 20000 samples with 
different configurations of defects for systems with spin 
concentrations p = 0.95 and p = 0.8, accordingly. 

We have analyzed the time dependence of the cumu- 
lant Uiify for samples with spin concentrations p = 0.8 
and clarified that in the time interval t G [10, 50] MCS/s 
the V^^t") is best fitted by power law with the dynamic 
exponent z = 2.068(24), corresponding to the pure Ising 
model [2J, |3J] and the influence of defects is developed 
for t > 400 MCS/s only. An analysis of the U2{t) slope 
measured in the interval t e [500, 950] MCS/s shows that 
the exponent d/z = 1.268(15) which gives z = 2.366(28). 
We have taken into account these dynamic crossover ef- 
fects for the analysis of the time dependence of magneti- 
zation and its derivative. So, the slope of magnetization 
measured in the interval t G [400, 950] MCS/s and its 
derivative over the interval t € [500, 950] MCS/s provides 
the exponents P/vz = 0.213(2) and l/i^z = 0.600(8) 
which give i^ = 0.704(18) and /3 = 0.365(8). The same 
analysis of the observable variables for samples with 
spin concentrations p — 0.95 leads to the value of ex- 
ponent d/z = 1.475(12), with z = 2.034(16), in the 
time interval t £ [10, 200] MCS/s and to the exponents 
d/z = 1.369(13), 13/iyz = 0.213(2), and 1/i^z = 0.600(8) 
in the time interval t e [550, 950] MCS/s, which give 
z = 2.191(21), 1/ = 0.704(18), and ^ = 0.365(8). 

For demonstration of crossover effects between the 
pure and the dilute regimes, we inserted in Fig. 2(a) and 
2(b) results of the magnetization measurements for pur e 
Ising model at the critical temperature Tc = 4.51142 [3a ]. 
Comparison of obtained curves in the insets confirms 
our conclusion that the influence of disorder on the non- 
equilibrium critical relaxation is developed for t > 550 
MCS/s for samples with p = 0.95 and for t > 400 MCS/s 
for samples with p — 0.8. 

In the next stage, we have considered the corrections to 
the scaling in order to obtain accurate values of the crit- 
ical exponents. We have applied the following expression 
for the observable X{t): 



Xit) = A,t'{l + B.J-^/') 



(16) 



where a; is a well-known exponent of corrections to scal- 
ing. Ax and B^ are fitting parameters, and an exponent 
S = —P/vz when X = m{t), S — d/z when X = U2{t), 
and S = 1/vz when X = c^t- lnm(f, r)|^^Q. This ex- 
pression reflects the scaling transformation in the criti- 
cal range of time-dependent corrections to scaling in the 
form of t^"^/^ to the usual form of corrections to scal- 
ing T^'^ in equilibrium state for time t comparable to 
the order parameter relaxation time t^ '^ £,^n,{k^) (ITJ . 
Field-theoretic estimate of the uj value gives w ~ 0.25(10) 
in the six-loop approximation [36|- Monte Carlo studies 
show that UJ ~ 0.370(63) from Ref. [13 and w ~ 0.26(13) 
fromRef. HI. 
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FIG. 2: Time evolution of the magnetization m{t) is plotted on a log-log scale for samples with spin concentrations p — 0.95 
(a) and p = 0.8 (b) at T = Tc{p) (curves 1) and at T = Tc{p) T AT with AT = 0.005 (curves 2 and 3). In insets curves 4 
correspond to pure Ising model. 
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FIG. 3: Time evolution of the logarithmic derivative of the magnetization 9r ln7n(t,r)|^^y with respect to t is plotted on a 
log-log scale for samples with spin concentrations p = 0.95 (a) and p = 0.8 (b). 
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FIG. 4: Time evolution of the cumulant U2{t) is plotted on a log-log scale at T = Tc{p) for samples with spin concentrations 
p = 0.95 (a) and p = 0.8 (b). 




0.225 0.240 



P/VZ " 





1.32 1.34 1 36 1 38 J/z ^ ^^ ^ ''2 



FIG. 5: Dependence of the mean square errors a of the fits for the magnetization (a), logarithmic derivative of the magnetization 
(b), and the cumulant (c) as a function of the exponents P/vz, l/vz, and d/ z for p — 0.8. 



TABLE I: Values of the exponents P/i/z, d/z, l/vz, and oj/z, 
corresponding minimal values of the mean-square errors for 
spin concentrations p = 0.95 and p = 0.80. 



p 


Exponent 


Mean 
value 


Approximation 
errors 


Statistical 
errors 


^/z 


0.95 


P/t.z 


0.244 


0.00011 


0.00131 


0.234 




d/z 


1.373 


0.00938 


0.00642 


0.092 




1/uz 


0.685 


0.00117 


0.00583 


0.181 


0.80 


Pjvz 


0.230 


0.00081 


0.00393 


0.275 




d/z 


1.359 


0.01209 


0.00785 


0.132 




\/vz 


0.661 


0.00418 


0.00700 


0.142 



We have used the least-squares method for the best 
approximation of the simulation data X{t) by the expres- 
sion in Eq. (J16p . Minimum of the mean square errors a 
of this fitting procedure determines the exponents 6 and 
io/ z. As example, for samples with spin concentrations 
p = 0.8, we plot in Fig. 5 the a for the magnetization (a) 
as a function of the exponent fi/vz for lo/ z — 0.275, loga- 
rithmic derivative of the magnetization (b) as a function 
of the exponent \/vz for w/z = 0.142, and the cumulant 
(c) as a function of the exponent d/z for uj/ z ~ 0.132. 
In Table I we present the computed values of the expo- 
nents ji/vz^ d/z, 1/j^z, and w/z, corresponding minimal 
values of the mean square errors a in these fits. The 
statistical errors for exponents are estimated by divid- 
ing all data into five data sets. On the base of these 
values of exponents and average value of w/z, we de- 
termine the final values of the critical exponents z — 
2.185(25), ^/v = 0.533(13), v = 0.668(14), ^ = 0.356(6), 
and a; = 0.369(96) for p = 0.95, and z = 2.208(32), 
^/v = 0.508(17), V = 0.685(21), /3 = 0.348(11), and 
a; = 0.404(110) for p = 0.8. 

The comparison of the obtained values of critical ex- 
ponents shows their belonging to the same class of uni- 
versal critical behavior of the diluted Ising model which 
can be characterized by the averaged critical exponents 
z = 2.196(17), V = 0.677(11), /3 = 0.352(5), and 
a; = 0.387(60). 



B. Evolution from a disordered state with mo <g 1 

In this part of the paper, we present the numerical 
investigations of the short-time critical dynamics of the 
3D site-diluted Ising model on the lattice with linear size 
L = 128, starting from a disordered state with small ini- 
tial magnetizations ttiq = 0.01, 0.02, and 0.03 for samples 
with spin concentration p — 0.8 only. For independent 
determination of the dynamic critical exponent z and 
the ratio of static exponents fi /v we investigate also a 
time dependence of the second moment of magnetization 



ra 



(2) 



(t) and the auto-correlation function A{t) for sys- 
tem, starting from a high-temperature initial state with 
Too = (in fact, with mo = 10~^). In accordance with 
Sec. II, a generalized dynamic scaling predicts in this case 
a power law evolution for the magnetization m(t), the 
second moment TO,"-'(t) and the autocorrelation function 
A(t) in the short-dynamic regime. 

Initial configurations for systems with the initial mag- 
netization Too were generated numerically. The initial 
magnetization has been prepared by flipping in an or- 
dered state a definite number of spins at randomly cho- 
sen sites in order to get the desired small value of too. 
Starting from those initial configurations, the system was 
updated with Metropolis algorithm at the critical tem- 
perature Tc — 3.49948(18), which has been determined 
in our paper [33] using Monte Carlo simulation of the 
3D site-diluted Ising model with p — 0.8 in equilibrium 
state. 

We measured the time evolution of the magnetization 
m(t) with values of the initial magnetization toq = 0.01, 
0.02, and 0.03, the second moment m^^)(i) and the au- 
tocorrelation function A{t) with mo = 0.0001 up to 
t = 1000 MCS/s. We show the obtained curves for 171(1) 
in Fig. 6, for m^^^t) in Fig. 7a, and for A{t) in Fig. 7b, 
which are plotted on a log-log scale. These curves were 
obtained by averaging over 4000 different samples with 
25 runs for each sample. We can see an initial increase 
of the magnetization, which is a very prominent phe- 
nomenon in the short-time critical dynamics. But in 
contrast to dynamics of the pure systems [2J|, we can 
observe the crossover from dynamics of the pure system 
on early times of the magnetization evolution from t ~ 15 




TABLE III: Values of the exponents for the 3D site-diluted 
Ising model to p = 0.80 obtained with the use of corrections 
to the scaling. 



t(MCS/s)' 



Exponent 


Value 


Lj/z 


6»'(mo = 0.03) 
9' {mo = 0.02) 
9'{mo = 0.01) 


0.104(12) 
0.117(10) 
0.118(10) 


0.074 
0.068 
0.096 


9' {mo -^ 0) 
C2(mo = 0) 
Ca{mo = 0) 


0.127(16) 
0.909(4) 
1.242(10) 


0.079 
0.112 
0.160 


(a;/z)av 

(^)av 


2.191(21) 
0.504(14) 
0.117(24) 
0.256(55) 





FIG. 6: Time evolution of the magnetization m,{t) for 
different values of the initial magnetization mo = 0.01 (1); 
0.02 (2); 0.03 (3), plotted on a log-log scale for samples with 
spin concentration p — 0.8. 



TABLE II: The initial slip exponent 9' measured by simula- 
tion of the 3D site-diluted Ising model with p = 0.80 for dif- 
ferent values of the initial magnetization m,o and exponents 
C2 and Ca for mo — 0. The value 9' {mo = 0) is the result of 
an extrapolation. 



mo 



0.03 

0.02 

0.01 





0.03 

0.02 

0.01 





C2 



l3/u 



t G [15, 601 



t e [5, 301 



0.1016(9) 
0.1031(10) 
0.1043(12) 
0.1057(17) 0.936(4) 



1.347(8) 2.065(14) 0.534(6) 



t e [300, 800] t g [150, 



0.083(3) 
0.099(5) 
0.105(9) 
0.122(11) 



0.859(5) 1.135(10) 2.387(20) 0.475(14) 



up to t ~ 60 MCS/s to dynamics of the disordered system 
with the influence of point-like defects in the time inter- 
val t € [300, 800] MCS/s. The same crossover phenomena 
were observed in evolution of the second moment of mag- 
netization ■m^^'>{t) and the autocorrelation function A{t). 
In the result of linear approximation of these curves in 
both the time intervals we obtained the values of the ex- 
ponent 6'{mQ) for initial states with toq = 0.01, 0.02, 
and 0.03 and the exponents C2 and Ca in accordance with 
relations in Eqs. ©, (0), and (HU]) (Table II). The final 
value of 0' is determined by extrapolation to rriQ — 0. 
Note that the similar crossover phenomena in the non- 
equilibrium critical relaxation of systems with quenched 
disorder have been revealed earlier in Ref. [31| by means 
of numerical simulations of the critical behavior of the 
3D Ising and XY models with linear defects. 

In the next stage, we have applied the procedure of 



corrections to the scaling determining by the expression 
([T6| for analysis of the observable m{t), m^-^^t), and 
A{t). We have used the least-squares method for the 
best approximation of the simulation data by the ex- 
pression (J16p . Minimum of the mean-square errors a of 
this fitting procedure determines the exponents 9' (mo), 
C2, and Ca with their respective tu/z. In Table III, we 
present the computed values of these exponents and the 
final value of 9' = 0.127(16) obtained by extrapolation 
of 9' {mo) for different values of the initial magnetization 
rriQ to mo = 0. In Table III, we also give the values of 
critical exponents z^ fijv-, and the average value of w and 
compare the values of these exponents to values of cor- 
responding exponents for the pure Ising model 241] . The 
obtained values agree quite well with results of simulation 
from an ordered state with toq ~ 1. 

Comparison of the value 9' = 0.127(16) to 9' = 0.10(2) 
from Ref. 27] also measured by simulation of the 3D site- 
diluted Ising model shows their not bad agreement within 
the limits of statistical errors of simulation and numerical 
approximations. In Ref. [27| the non-equilibrium relax- 
ation of the magnetization m(i) has been investigated 
from the initial random spin configurations with mean 
magnetization niQ — 0.01 only for samples with different 
spin concentrations p = 0.499, 0.6, 0.65, and 0.8 and with 
linear sizes L = 8, 16, 32 and 64. However, in accordance 
with Ref. [2J], the initial slip exponent 9' must be deter- 
mined in the asymptotical limit with rtiQ ^ on the basis 
of results of 'm{t) computing for a few small values of the 
initial magnetization mo. Furthermore, as it follows from 
Eq. (5), the magnetization undergoes a powerlaw initial 
increase characterized by 9' for sufficiently small t^°^^mo. 
For rriQ and t not too small, the powerlaw behavior will 
be modified. For strongly diluted systems which are also 
considered in Ref. [23|, the influence of disorder is ob- 
served for longer times than for weakly diluted systems. 
Therefore, the use of the same value of the initial mag- 
netization for determination of the initial slip exponent 
9' for both weakly and strongly diluted systems is un- 
justified. It should be noted that in Ref. ^3 the data 
analysis of m{t) for samples with different spin concen- 
trations p was carried out with the use of the corrections 
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FIG. 7: Time evolution of the second moment m,"' (a) and the correlation function A{t) (b) for L — 128 with the initial 
magnetization mo = 0.0001. 



to scaling procedure during realization of which the uni- 
versal value of the dynamic critical exponent z = 2.62(7) 
obtained in Ref. [3^ was applied. But this value z is 
inconsistent with values calculated both in the present 
paper and in Ref. [39| in the three-loop approximation 
of the field-theoretic RG description, and with experi- 
mentally measured value z — 2.18(10) for weakly diluted 
Ising magnet Feo.9Zno.1F2 from Ref. j|40]. 

Obtained in present paper is the asymptotic value 
6'{mQ ~^ 0) — 0.127(16) that demonstrates that it is 
larger than 6'{mo = 0.01) from Ref. i27i]. It is explained 

by the revealed tendency that d'irriQ ) > O'irriQ ) for 
the initial magnetization which are in the following cor- 
respondence with each other as ttIq < rri n ■ Therefore, 
the distinguished good agreement in Ref. |27| of the ob- 
tained value e' = 0.10(2) with 9' ~ 0.0867 from Ref. ^ 
is unjustified. The results of investigations carried out 
in this paper give reasons to consider that the value of 
the initial slip exponent 9' = 0.127(16) is more realistic 
for description of non-equilibrium critical relaxation of 
the 3D weakly diluted Ising-like systems which is larger 
than the value of exponent 6' = 0.108(2) for the pure 3D 
Ising systems [1^ |24| rather than smaller as predicted by 
results from Refs. [23| and |29l |. 

We have realized a field-theoretic renormalization- 
group description of non-equilibrium critical relaxation 
for directly three-dimensional diluted Ising model and 
calculated the initial slip exponent 9' in two-loop ap- 
proximation without using the e-expansion method. As 
a result, we have obtained 



9' = — 5* + 0.125W* - 0.123968(.g*)2 + 
6 

+ 0.14680608g*i;* - 0.0156245(l;*)^ 



(17) 



where g* and v* are values of the vertexes describ- 
ing interaction of the order parameter fluctuations in 
the fixed-point of the renormalization-group equations 
I41I l42i. For further calculations, we use the FP with 



g* = 2.2514(42), v* = -0.7049(13) which determines the 
critical behavior of the 3D dilute Ising model. The co- 
ordinate of this FP has been obtained in our paper [33 
as average of numerical values g*, v* which were calcu- 
lated with the use of different methods of resummation 
technique. It is well known that the series expansions for 
the critical exponents exhibit factorial divergence, but 
they can be considered in an asymptotic context. In or- 
der to obtain physically reasonable values of the critical 
exponents for 3D systems, special methods for the sum- 
mation of asymptotic series have been developed f3y, |43| - 
47], the most effective being the Pade-Borel, Pade-Borel- 
Leroy (PBL), and conformal mapping techniques. We 
employ the PBL resummation method extended to the 
two-parameter case. The PBL method is a generaliza- 
tion of the Pade-Borel method with the integral Borel 
transformation 
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00 
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(18) 



where the value of parameter b = 2.221426 was chosen in 
Ref. [331 from convergence analysis of the test series for 
exactly solvable problem of calculating the anharmonic 
oscillator energy with the asymptotic convergence of the 
series, which is similar to the series for the RG /? and 7 
functions in the theory of critical phenomena. As a result, 
the following value 9' — 0.1203 has been calculated which 
very well agrees with our results of simulation. 



IV. ANALYSIS OF RESULTS AND 
CONCLUSIONS 

In a summary Table IV, we present the values of crit- 
ical exponents z, 0', /3/i^, u^ /3, and w obtained in this 



TABLE IV: Values of the obtained critical exponents and comparison to other results of Monte Carlo simulations (MC), 
field-theoretical method with fixed-dimension d = 3 expansion (FTM), and experimental (EXP) investigations 





z 


e' 


P/v 


V 


/? 


Ix) 


p = 0.95, mo = 1 




2.185(25) 




0.533(13) 


0.668(14) 


0.356(6) 


0.369(96) 


p = 0.80, mo = 1 




2.208(32) 




0.508(17) 


0.685(21) 


0.348(11) 


0.404(110) 


p = 0.80, mo < 1 




2.191(21) 


0.127(16) 


0.504(14) 






0.256(55) 


Pelissetto, Vicari, 2000, (Ref. [36]) 


(FTM) 






0.515(15) 


0.678(10) 


0.349(5) 


0.25(10) 


Prudnikov, et al, 2006, (Ref. [39]) 


(FTM) 


2.1792(13) 












Rosov, et al, 1988, FepZni-pF2 p = 0.9, (Ref. [56]) 


(EXP) 










0.350(9) 




Rosov, et al, 1992, Fe^Zni-pFa p = 0.9, (Ref. [40]) 


(EXP) 


2.18(10) 












Slanic, et al, 1999, FepZni_pF2 p = 0.93, (Ref. [54]) 


(EXP) 








0.70(2) 






Prudnikov, Vakilov, 1992, p = 0.95, 




2.19(7) 












p = 0.80, 




2.20(8) 












p = 0.60, 




2.58(9) 












p = 0.40, (Ref. [52]) 


(MC) 


2.65(12) 












Heuer, 1993, p = 0.95 




2.16(1) 




0.49(2) 


0.64(2) 


0.31(2) 




p = 0.90 




2.232(4) 




0.48(2) 


0.65(2) 


0.31(2) 




p = 0.80, 




2.38(1) 




0.51(2) 


0.68(2) 


0.35(2) 




p = 0.60, (Refs. [50], [51]) 


(MC) 


2.93(3) 




0.45(2) 


0.72(2) 


0.33(2) 




Wiseman, Domany, 1998, p = 0.80, 








0.505(2) 


0.682(2) 






p = 0.60, (Ref. [48]) 


(MC) 






0.437(21) 


0.717(6) 






Ballesteros, et al, 1998, p = 0.90 -^ 0.40 (Ref. [37]) 


(MC) 






0.519(8) 


0.684(5) 


0.355(3) 


0.370(63) 


Parisi, ei a/., 1999, p = 0.90 -;- 0.40 (Ref. [38]) 


(MC) 


2.62(7) 










0.50(13) 


Calabrese, et al, 2003, p = 0.80 (Ref. [49]) 


(MC) 






0.518(5) 


0.683(3) 


0.354(2) 




Murtazaev, et al, 2004, p = 0.95, 










0.646(2) 


0.306(3) 




p = 0.9. 










0.664(3) 


0.308(3) 




p = 0.8, 










0.683(4) 


0.310(3) 




p = 0.6, (Ref. [53]) 


(MC) 








0.725(6) 


0.349(4) 




Schehr, Paul, 2006, (Ref. [27]) 


(MC) 




0.10(2) 










Hasenbusch, et al, 2007, p = 0.8 (Ref. [55]) 


(MC) 


2.35(2) 












Prudnikov, et a/., 2007, p = 0.95 -^ 0.80, 








0.532(12) 


0.693(5) 




0.26(13) 


p = 0.60 H- 0.50, (Ref. [33]) 


(MC) 






0.524(13) 


0.731(11) 




0.28(15) 



paper by comprehensive Monte Carlo simulations of the 
short-time critical evolution of the 3D site-diluted Ising 
model both from an ordered initial state with mo — 1 for 
samples with spin concentrations p = 0.95 and 0.8 and 
from a disordered initial states with ttiq ^ 1 for samples 
with spin concentration p ~ 0.8. For comparison, we give 
in Table IV the results of calculation of these exponents 
by the field-theoretical method with fixed-dimension d = 
3 expansion [3a, [33, the results of experimental inves- 
tigations of the Ising-like magnets 40l 54 561 . and the 
results of numerical studies [23, |3J, [S^, [33, l48l - [53l ISSf . 
As shown in Table IV, our values of exponents are in 
good agreement within the limits of statistical errors of 
simulation and numerical approximations with results of 
the field-theoretical description of the statics in six-loop 
approximation 13 61 and critical dynamics in three-loop 
approximation [39[ and with the results of experimen- 
tal investigations of the static [54;, ^6] and dynamic fioj 
critical behaviors of weakly diluted Ising-like magnets. 

Comparison to results of Monte Carlo simulations 
shows that our static exponents f^jv, v, and ^ agree 
well with those exponents measured in equilibrium for 
weakly diluted systems in the most cited pap ers and for 
systems with wide dilution range in Ref. [37| where the 
critical exponents were obtained for p < 0.8 as dilution- 



independent after a proper infinite volume extrapola- 
tion with taking into account the leading corrections-to- 
scaling terms. However, it was found that the case with 
•p — 0.9 falls out from this dilution-independent scheme 
of fits with common exponent a; — 0.37 for samples with 
different spin concentrations. Authors draw a conclusion 
that the p = 0.9 data seem to be still crossing over from 
the pure Ising fixed point to the diluted one. Most of 
the computations have been carried out at p = 0.8 as in 
this case the scaling corrections are very small and the 
results even in small lattices are stable. So, in Ref. j49{ , 
it was shown that for case with p = 0.8, the observed 
corrections to scaling could be the next-to-leading with 
a;2 "^ 0.80. However, our present investigation and the 
results of computations in equilibrium [33] show that the 
systems with p = 0.95 and p = 0.8 are characterized by 
close agreement of the critical exponent values and they 
belong to the same class of universal critical behavior of 
the site-diluted Ising model with the averaged critical ex- 
ponents V = 0.677(11), /3 = 0.352(5), and w = 0.387(60). 

Now we compare the values of the dynamic critical ex- 
ponent z obtained in this paper by the short-time dynam- 
ics method to the results of Monte Carlo simulations of 
the critical dynamics in equilibrium, realized in Refs. |5l| 
and [5 51 , to the results of non-equilibrium studies of the 
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susceptibility in Refs. [Sg] and 55 1, and to the results 
of Monte Carlo renormalization group application to de- 
scription of the site-diluted Ising model relaxation from 
the ordered initial state with toq = 1 in Ref. 52] . Values 
of z from paper [51] agree rather well with our results 
only for weakly diluted systems with p > 0.9, while a 
noticeable difference between the results is observed for 
strongly disordered systems. Starting from the universal- 
ity concept for critical behavior of diluted Ising systems 
and that the asymptotic value of z is independent of the 
degree of dilution, the author in Ref. [5l| obtained the 
asymptotic value z — 2.4(1) using the effective values 
of the exponent listed in Table IV. The off-equilibrium 
critical dynamics of the 3D Ising model with the spin 
concentration varying in a wide range was analyzed in 
Ref. [38.]. Also assuming that the critical behavior of 
diluted Ising systems is universal under dilution, the au- 
thors obtained the asymptotic value oiz = 2.62(7) taking 
into account the leading corrections to the scaling depen- 
dence for the dynamical susceptibility. In this case, the 
value of the exponent w = 0.50(13) obtained in Ref. [S^ 
is strongly inconsistent with w — 0.25(10) from the field 
theory calculations |36| and not so well agreement with 
w = 0.37(6) from Monte Carlo results in Ref. [33|. In 
the approximations realized in Ref. jSg], the results for 
weakly diluted systems were characterized by the largest 
errors. In Ref. [55j], it was carried out the Metropolis dy- 
namics in equilibrium for site-diluted Ising model with 
p = 0.8, 0.85, and 0.65. For case with p = 0.8, authors 
investigated in detail the scaling corrections which were 
the next-to-leading with 102 = 0.82(8) and gave the expo- 
nent z = 2.35(2). Investigations for other values of p did 
not permit to determine z in these systems accurately. 
Also, in Ref. [5a], investigated was the off-equilibrium 
relaxational critical dynamics in the site-diluted Ising 
model dX p — 0.8. The results show that equilibrium 
estimate z = 2.35(2) is perfectly consistent with the 
off-equilibrium MC data. Authors did not observe a 
large-time scaling corrections proportional to t~^l^\ in- 
stead, their data show corrections that are proportional 
to t~^'^l^ with the static correction-to-scaling exponents 
w = 0.29(2) and uj2 = 0.82(8). We have some doubts 
in validity of z = 2.35(2). As it was shown in [33, the 
realization of correction to scaling procedure demands 
the six simulation data points at least for their approx- 
imation by four-parameter function such as in Eg. (J16p 
for lattices with L > imin- Whereas in paper [53|, the 
asymptotical value z = 2.35(2) was obtained with the 
use of four or five data points as it was demonstrated in 
Figs. 1-4 in [53 • It is necessary to use additional data 
points for lattices with i > 64. 

The early results of our numerical investigations of the 
critical dynamics for diluted Ising systems in Ref. 52] by 
Monte Carlo renormalization-group method show very 
good agreement with our present results for weakly di- 
luted systems, while a noticeable difference between the 
results is observed for strongly diluted systems. We are 
planning to continue the Monte Carlo study of critical 



behavior of the site-diluted Ising model by short-time 
dynamics method with p — 0.6 and 0.5 focusing on the 
problem of dilution independence of asymptotic charac- 
teristics. 

The present results of Monte Carlo investigations allow 
us to recognize that the short-time dynamics method is 
reliable for the study of the critical behavior of the sys- 
tems with quenched disorder and is the alternative to 
traditional Monte Carlo methods. But in contrast to 
studies of the critical behavior of the pure systems by 
the short-time dynamics method [23, [2J| , in case of the 
systems with quenched point-like disorder after the mi- 
croscopic time tmic — 5 -^ 10 MCS/s, there exist three 
stages of dynamic evolution. For systems starting from 
the ordered initial states (ttiq = 1) in the time interval 
of 10-50 MCS/s, the power-law dependence is observed 
at the critical point for Binder cumulant U2{t), which 
is similar to that in the pure system. In the time in- 
terval [400,950], the power- law dependences are observed 
in the critical point for the magnetization 'm{t), the log- 
arithmic derivative of the magnetization, and cumulant 
U2{t) which are determined by the influence of disorder. 
However, careful analysis of the slopes for m{t) and U2{t) 
reveals that a correction to scaling should be considered 
in order to obtain accurate results. The dynamic and 
static critical exponents were computed with the use of 
the corrections to scaling, which demonstrate their good 
agreement with results of the field-theoretic description 
of the critical behavior of these models with disorder. 
In the intermediate time interval of 100-400 MCS/s, the 
dynamic crossover behavior is observed from the critical 
behavior typical for the pure systems to behavior deter- 
mined by the influence of disorder. 

The investigation of the critical behavior of the Ising 
model with defects starting from the disordered initial 
states with mp <C 1 also has revealed three stages of 
the dynamic evolution. It was shown that the power-law 
dependences for the magnetization m(t), the second mo- 
ment m!''^\t), and the autocorrelation A{t) are observed 
in the critical point, which are typical for the pure system 
in the common interval [5,60] for observable quantities 
and for the disordered system in the common interval 
[150,800]. In the intermediate time interval the crossover 
behavior is observed in the dynamic evolution of the sys- 
tem. The obtained values of exponents demonstrate good 
agreement within the limits of statistical errors of simula- 
tion and numerical approximations with results of simu- 
lation of the pure Ising model by the short-time dynamics 
method for the first time interval and with our results of 
simulation of the critical relaxation of this model from 
the ordered initial state. 
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